Module MCMC_subroutines

  ! Module contains pdf subroutine, randomn function and random_seed subroutine, inverse and (log sqrt) determinant fuctions, Kronecker product (not used in this version), gamma function, bessel
  use MCMC_utils, only        : pdf_calc, random_normal, init_random_seed, &
       Ainv, det, kron, gamma, rkbesl

  Implicit none

Contains

  subroutine setup
    use MCMC_data

    integer :: ii, jj, xi, yi

    ! Set first value of each parameter from prior values
    x=x_ap
    pdf_param1_current = pdf_param1
    pdf_param2_current = pdf_param2
    sigma_y_current = sigma_y
    sigma_ys_current = sigma_ys
    tau_current = tau
    nu_current = nu
    rho_current = rho

    ! Create autoregressive AR(1) temporal autocorrelation matrix

    do ii=1,dim2
       y_error_t(T_indices(:,ii)) = sigma_y_current(ii)
    enddo

    !$OMP PARALLEL DO default(none), &
    !$OMP             private(autocorr_vec), &
    !$OMP             shared(T_current, y_error_t, deltatime, tau_current, nmeasuretotal)
    do ii=1,nmeasuretotal
       autocorr_vec = y_error_t(ii)*y_error_t*exp(-1*(deltatime(:,ii)/tau_current))
       T_current(ii,:) = autocorr_vec
       T_current(:,ii) = autocorr_vec
    enddo
    !$OMP END PARALLEL DO

    ! Compute inverse of T

    Tinv_current = Ainv(T_current)

    ! find number of elements corresponding to each column of T_indices in order to scale determinant

    do ii = 1,dim2
       num_T_indices(ii) =maxloc(abs(T_indices(:,ii)), 1)
    enddo

    ! Create Matern covariance for spatial correlation

    nu_double = nu_current
    rho_double = rho_current
    distance_double = distance
    nb = int(nu_double)
    alpha = nu_double - nb
    nb = nb+1
    ize = 2

    arg = sqrt(2*nu_double)*distance_double/rho_double

    call gamma(nu_double, ga_double)
    ga = ga_double

    !$OMP PARALLEL DO default(none), &
    !$OMP             collapse(2), &
    !$OMP             private(arg_dum, k_arg_dum), &
    !$OMP             shared(S_current, numsites, arg, sigma_ys_current), &
    !$OMP             shared(ncalc,ize,nb,alpha,nu_current,distance,rho_current,ga)
    do ii = 1,numsites
       do jj = 1,numsites
          arg_dum = arg(ii,jj)
          if (arg_dum==0) then
             S_current(ii,jj) = 1*sigma_ys_current(ii)**2
          else
             call rkbesl(arg_dum, alpha, nb, ize, k_arg_dum, ncalc)
             k_arg_dum = k_arg_dum*exp(-1*arg(ii,jj))
             S_current(ii,jj) = sigma_ys_current(ii)*sigma_ys_current(jj)*(1/(2**(nu_current-1)*ga)* &
                  (sqrt(2*nu_current)*distance(ii,jj)/rho_current)**nu_current*k_arg_dum)
          endif
       enddo
    enddo
    !$OMP END PARALLEL DO

    ! Compute inverse of S
    Sinv_current = Ainv(S_current)

    ! Calulate R using the Kronecker product
    if (kron_flag==0) then
       call kron(Rinv_current, Sinv_current, Tinv_current)
    elseif (kron_flag==1) then
       call kron(Rinv_current, Tinv_current, Sinv_current)
    endif

    ! Constants for Intel fortran matrix multiplication routines
    aa = 1
    bb = 0

    ! Compute (log sqrt) determinant of R using properties of the Kronecker product
    detval_S_current = det(S_current)*nmeasuretotal
    detval_T_current = det(T_current)*numsites
    detval_current = detval_S_current + detval_T_current

    ! Calclulate the inverse of the measurement uncertainty nugget term (assumed diagonal)
    Dinv = Ainv(D)

    ! Intial value of model simulated mole fractions
    y_current = matmul(H,x)

    ! Subsample at measurement times
    y_small = y_current(timeindex_nonzero)

    ! Initial value of the difference between observations and model mole fractions at observation times
    m0 = z - y_small

    ! Initialise m1 and derived values (used in compute_y)
    m1 = m0
    do ii=1,nmeasure
       dum2(ii) = m1(ii)*Dinv(ii,ii)*m1(ii)
    end do
    m1T = sum(dum2)

    ! Set up timeindex array

    timeindex(:) = 0
    do ii=1,nmeasure
       timeindex(timeindex_nonzero(ii)) = ii
    end do

    ! Compute (z-y)^T*Dinv*(z-y)
    call sgemm ('N', 'N', nmeasure, 1, nmeasure, aa, Dinv, nmeasure, m0, nmeasure, bb, C2, nmeasure)

    m0T = dot_product(m0,C2)

    ! Compute difference between sampled value of model mole fractions and model derived ones based on value of x
    n0 = matmul(H,x) - y_current
    n1 = n0

    ! Compute (Hx-y)^T*Rinv*(Hx-y) by reshaping into an array and using properties of the Kronecker product
    C = reshape(matmul(matmul(Tinv_current,reshape(n0, (/nmeasuretotal, &
                numsites/))),Sinv_current),(/nmeasuremax/))

    n0T = dot_product(n0,C)


    ! Compute P0 for starting state
    do xi=1,statesize
       call pdf_calc(x(xi), pdf_param1(xi), pdf_param2(xi), x_pdf(xi), p1)
       p0_x(xi) = p1

       call pdf_calc(pdf_param1_current(xi), pdf_param1(xi), pdf_hyperparam1(xi), pdf_param1_pdf(xi), p1)
       p0_pdf_param1(xi) = p1

       call pdf_calc(pdf_param2_current(xi), pdf_param2(xi), pdf_hyperparam2(xi), pdf_param2_pdf(xi), p1)
       p0_pdf_param2(xi) = p1
    enddo

    do yi = 1,dim2
       call pdf_calc(sigma_y_current(yi), sigma_y(yi), sigma_y_hyperparam(yi), sigma_y_pdf(yi), p1)
       p0_sigma_y(yi) = p1
    enddo

    do yi = 1,numsites
       call pdf_calc(sigma_ys_current(yi), sigma_ys(yi), sigma_ys_hyperparam(yi), sigma_ys_pdf(yi), p1)
       p0_sigma_ys(yi) = p1
    enddo

    call pdf_calc(tau_current, tau_hyperparam1, tau_hyperparam2, tau_pdf, p1)
    p0_tau = p1

    call pdf_calc(nu_current, nu_hyperparam1, nu_hyperparam2, nu_pdf, p1)
    p0_nu = p1

    call pdf_calc(rho_current, rho_hyperparam1, rho_hyperparam2, rho_pdf, p1)
    p0_rho = p1

    ! Initialize arrays that will store values and acceptance ratios
    x_it(1,:)=x
    pdf_param1_it(1,:) = pdf_param1
    pdf_param2_it(1,:) = pdf_param2
    sigma_y_it(1,:) = sigma_y
    sigma_ys_it(1,:) = sigma_ys
    tau_it(1,1) = tau
    nu_it(1,1) = nu
    rho_it(1,1) = rho
    y_it(1,:) = y

    reject=0
    reject_vector(:)=0
    reject_sigma_y=0
    reject_sigma_y_vector(:)=0
    reject_sigma_ys = 0
    reject_sigma_ys_vector(:)=0
    reject_tau=0
    reject_nu=0
    reject_y = 0
    reject_y_vector(:) = 0

    call init_random_seed(fixed_seed=1234) ! For repeatible testing, use a fixed seed every time we run
    !call init_random_seed()              ! Ensure random number generation starts from new point each time program is run
    ! Random seed only needs to be called once in a program.
    ! It is used to generate the first random number.
    ! Any random numbers generated after the first will use the previous random number as its seed.

  end subroutine setup

  subroutine compute_x(stepsize_pdf_param1,pdf_param1,pdf_param1_pdf,&
    pdf_param1_current,stepsize_pdf_param2,pdf_param2,pdf_param2_pdf,&
    pdf_param2_current,stepsize,H,pdf_hyperparam1,dy,&
    pdf_hyperparam2,x,x_pdf,n1,n0,nmeasuretotal,nmeasuremax,&
    numsites,aa,Tinv_current,dum3,Sinv_current,bb,C,p0_x,&
    p0_pdf_param1,p0_pdf_param2,n0T,reject,reject_vector,burn_in_done)

    ! input only
    integer, intent(in) :: nmeasuretotal,numsites,nmeasuremax,pdf_param1_pdf,&
      pdf_param2_pdf,x_pdf
    real, intent(in) :: stepsize_pdf_param1,stepsize_pdf_param2,stepsize,pdf_hyperparam1,&
      pdf_hyperparam2,pdf_param1,pdf_param2,aa,bb
    real, dimension(:,:), intent(in) :: Tinv_current, Sinv_current
    real, dimension(:), intent(in) :: H
    logical, intent(in) :: burn_in_done

    ! inout
    real, intent(inout) :: pdf_param1_current, pdf_param2_current, x, p0_x,&
      p0_pdf_param1, p0_pdf_param2,n0T
    real, dimension(:), intent(inout) :: n0
    integer, intent(inout) :: reject, reject_vector

    ! really locals (temporary values)
    real, dimension(:), intent(inout) :: C,dy,n1
    real, dimension(:,:), intent(inout) :: dum3

    real :: dpdf_param1, dpdf_param2, dx, pdf_param1_new, pdf_param2_new, p1_pdf_param1,&
      p1_pdf_param2,p1_x,n1T,pT,randomu

    ! Calculate new values of x, pdf_param1, pdf_param2 using random number generator
    dpdf_param1 = random_normal()
    dpdf_param1 = dpdf_param1*stepsize_pdf_param1
    pdf_param1_new = pdf_param1_current + dpdf_param1

    dpdf_param2 = random_normal()
    dpdf_param2 = dpdf_param2*stepsize_pdf_param2
    pdf_param2_new = pdf_param2_current + dpdf_param2

    dx = random_normal()
    dx = dx*stepsize
    dy=H(:)*dx

    ! Compute P1 for each parameter
    call pdf_calc(pdf_param1_new, pdf_param1, pdf_hyperparam1, pdf_param1_pdf, p1_pdf_param1)
    call pdf_calc(pdf_param2_new, pdf_param2, pdf_hyperparam2, pdf_param2_pdf, p1_pdf_param2)

    call pdf_calc((x+dx), pdf_param1_new, pdf_param2_new, x_pdf, p1_x)  

    ! Compute n1T = (y-Hx)^T*Rinv*(y-HX) for new value of x
    n1=(n0+dy)

    call sgemm ('N', 'N', nmeasuretotal, numsites, nmeasuretotal, aa, Tinv_current, nmeasuretotal, &
         n1, nmeasuretotal, bb, dum3, nmeasuretotal)
    call sgemm ('N', 'N', nmeasuretotal, numsites, numsites, aa, dum3, nmeasuretotal, Sinv_current, &
         numsites, bb, C, nmeasuretotal)

    n1T = dot_product(n1,C)

    ! Compute P1/P0
    pT=alog(p1_x*p1_pdf_param1*p1_pdf_param2/(p0_x*p0_pdf_param1*p0_pdf_param2))-0.5*(n1T - n0T)

    call random_number(randomu)              ! Generates uniformly distributed random number

    if(alog(randomu) .le. pT) then
       !;ACCEPT
       p0_x=p1_x
       p0_pdf_param1 = p1_pdf_param1
       p0_pdf_param2 = p1_pdf_param2

       x=x + dx
       pdf_param1_current = pdf_param1_new
       pdf_param2_current = pdf_param2_new

       n0=n1
       n0T=n1T

    else
       !;REJECT
       if(burn_in_done) then
          reject=reject + 1
          reject_vector=reject_vector + 1
       endif
    endif

  end subroutine compute_x

  subroutine record_all(It)
    use MCMC_data

    integer, intent(in) :: It

    real :: t1

    if(mod(It+1,100) .eq. 0) then           ! Print acceptance ratio for every 1000th iteration to see if everything looks ok
       call system_clock(count_1,count_rate,count_max)
       t1=real(count_1)/real(count_rate)
       write(*,*) "Time: ",t1-t0, "Iters:", It - burn_in
       write(*,*) 'x acceptance:        ', 1. - real(reject)/real(It-burn_in)/real(statesize)
       write(*,*) 'sigma_y acceptance:  ', 1. - real(reject_sigma_y)/real(It-burn_in)/real(dim2)

       write(*,*) 'sigma_ys acceptance: ', 1. - real(reject_sigma_ys)/real(It-burn_in)/real(numsites)

       write(*,*) 'tau acceptance:      ', 1. - real(reject_tau)/real(It-burn_in)
       write(*,*) 'nu_rho acceptance:   ', 1. - real(reject_nu)/real(It-burn_in)
       write(*,*) 'y acceptance:        ', 1. - real(reject_y)/real(It-burn_in)/real(nmeasuremax)
    endif
    x_it(it-burn_in,:)=x                    ! Record x value for It > burn_in
    pdf_param1_it(it-burn_in,:) = pdf_param1_current
    pdf_param2_it(it-burn_in,:) = pdf_param2_current
    sigma_y_it(it-burn_in,:)=sigma_y_current   ! Record sigma_y value for It > burn_in
    sigma_ys_it(it-burn_in,:)=sigma_ys_current   ! Record sigma_ys value for It > burn_in
    tau_it(it-burn_in,:)=tau_current        ! Record tau value for It > burn_in
    nu_it(it-burn_in,:)=nu_current        ! Record nu and rho value for It > burn_in
    rho_it(it-burn_in,:)=rho_current
    y_it(it-burn_in,:)=y_current   ! Record y value for It > burn_in

  end subroutine record_all

  subroutine compute_sigma_y(sigma_y_current,sigma_y,sigma_y_hyperparam,&
    sigma_y_pdf,stepsize_sigma_y,Tinv_new,Tinv_current,T_indices,nmeasuretotal,&
    numsites,nmeasuremax,aa,bb,dum3,Sinv_current,num_T_indices,&
    detval_T_current,detval_S_current,p0_sigma_y,n0T,n0,C,&
    detval_current,reject_sigma_y,reject_sigma_y_vector,burn_in_done)

    ! input only
    integer, intent(in) :: sigma_y_pdf,nmeasuretotal,&
      numsites,nmeasuremax,num_T_indices
    real, intent(in) :: sigma_y,sigma_y_hyperparam,aa,bb,detval_S_current,&
      stepsize_sigma_y
    real, dimension(:,:), intent(in) :: Sinv_current
    real, dimension(:), intent(in) :: n0
    logical, intent(in) :: burn_in_done
    integer, dimension(:), intent(in) :: T_indices

    ! inout
    real, intent(inout) :: sigma_y_current,detval_T_current,p0_sigma_y,n0T,&
      detval_current
    real, dimension(:,:), intent(inout) :: Tinv_current
    real, dimension(:), intent(inout) :: C
    integer, intent(inout) :: reject_sigma_y, reject_sigma_y_vector

    ! really locals (temporary values)
    real, dimension(:,:), intent(inout) :: Tinv_new, dum3

    real :: sigma_y_new,dsigma_y, p1_sigma_y, n1T,detval_T_new,detval_new,&
      randomu,pT

    ! Generate new value of sigma_y
    sigma_y_new = sigma_y_current + random_normal()*stepsize_sigma_y

    ! Calculate P1 for new value of sigma_y
    call pdf_calc(sigma_y_new, sigma_y, sigma_y_hyperparam, sigma_y_pdf,&
      p1_sigma_y)

    ! Calculate new value of inverse by scaling old value
    Tinv_new = Tinv_current
    Tinv_new(T_indices(:),:) = 1/(sigma_y_new/sigma_y_current)*Tinv_current(T_indices(:),:)
    Tinv_new(:,T_indices(:)) = 1/(sigma_y_new/sigma_y_current)*Tinv_new(:,T_indices(:))

    ! calculate new value of n1T = (y-Hx)^T * Rinv * (y-Hx)

    call sgemm ('N', 'N', nmeasuretotal, numsites, nmeasuretotal, aa, Tinv_new, nmeasuretotal, n0, &
         nmeasuretotal, bb, dum3, nmeasuretotal)
    call sgemm ('N', 'N', nmeasuretotal, numsites, numsites, aa, dum3, nmeasuretotal, Sinv_current, &
         numsites, bb, C, nmeasuretotal)

    n1T = dot_product(n0,C)

    ! caluclate new determinant by scaling
    detval_T_new = (num_T_indices*alog(sigma_y_new/sigma_y_current))*numsites + detval_T_current

    detval_new = detval_S_current + detval_T_new

    ! Compute P1/P0
    pT=alog(p1_sigma_y/p0_sigma_y) - detval_new + detval_current - 0.5*(n1T - n0T)

    call random_number(randomu)              ! Generates uniformly distributed random number

    if(alog(randomu) .le. pT) then
       !;ACCEPT
       p0_sigma_y=p1_sigma_y
       sigma_y_current = sigma_y_new
       Tinv_current = Tinv_new
       detval_current = detval_new
       detval_T_current = detval_T_new
       n0T=n1T
    else
       !;REJECT
       if(burn_in_done) then
         reject_sigma_y=reject_sigma_y+ 1
         reject_sigma_y_vector=reject_sigma_y_vector+ 1
       endif
    endif

  end subroutine compute_sigma_y

  subroutine compute_sigma_ys(yi,sigma_ys_current,sigma_ys,sigma_ys_hyperparam,&
    sigma_ys_pdf,stepsize_sigma_ys,Sinv_new,Sinv_current,nmeasuretotal,&
    numsites,nmeasuremax,aa,bb,dum3,Tinv_current,&
    detval_S_current,detval_T_current,p0_sigma_ys,n0T,n0,C,&
    detval_current,reject_sigma_ys,reject_sigma_ys_vector,burn_in_done)

    ! input only
    integer, intent(in) :: yi,sigma_ys_pdf,nmeasuretotal,&
      numsites,nmeasuremax
    real, intent(in) :: sigma_ys,sigma_ys_hyperparam,aa,bb,detval_T_current,&
      stepsize_sigma_ys
    real, dimension(:,:), intent(in) :: Tinv_current
    real, dimension(:), intent(in) :: n0
    logical, intent(in) :: burn_in_done

    ! inout
    real, intent(inout) :: sigma_ys_current,detval_S_current,p0_sigma_ys,n0T,&
      detval_current
    real, dimension(:,:), intent(inout) :: Sinv_current
    real, dimension(:), intent(inout) :: C
    integer, intent(inout) :: reject_sigma_ys,reject_sigma_ys_vector

    ! really locals (temporary values)
    real, dimension(:,:), intent(inout) :: Sinv_new, dum3

    real :: sigma_ys_new,dsigma_ys, p1_sigma_ys, n1T,detval_S_new,detval_new,&
      randomu,pT

    ! Generate new value of sigma_ys
    sigma_ys_new = sigma_ys_current + random_normal()*stepsize_sigma_ys

    ! Calculate P1 for new value of sigma_ys
    call pdf_calc(sigma_ys_new, sigma_ys, sigma_ys_hyperparam, sigma_ys_pdf,&
      p1_sigma_ys)

    ! Calculate new value of inverse by scaling old value
    Sinv_new = Sinv_current
    Sinv_new(yi,:) = 1/(sigma_ys_new/sigma_ys_current)*Sinv_current(yi,:)
    Sinv_new(:,yi) = 1/(sigma_ys_new/sigma_ys_current)*Sinv_new(:,yi)

    ! calculate new value of n1T = (y-Hx)^T * Rinv * (y-Hx)

    call sgemm ('N', 'N', nmeasuretotal, numsites, nmeasuretotal, aa,&
      Tinv_current, nmeasuretotal, n0, nmeasuretotal, bb, dum3, &
      nmeasuretotal)
    call sgemm ('N', 'N', nmeasuretotal, numsites, numsites, aa, dum3, &
      nmeasuretotal, Sinv_new, numsites, bb, C, nmeasuretotal)

    n1T = dot_product(n0,C)

    ! caluclate new determinant by scaling
    detval_S_new = alog(sigma_ys_new/sigma_ys_current)*nmeasuretotal + &
      detval_S_current

    detval_new = detval_T_current + detval_S_new

    ! Compute P1/P0
    pT=alog(p1_sigma_ys/p0_sigma_ys) - detval_new + detval_current - &
      0.5*(n1T - n0T)

    call random_number(randomu) ! Generates uniformly distributed random number

    if(alog(randomu) .le. pT) then
       !;ACCEPT
       p0_sigma_ys=p1_sigma_ys
       sigma_ys_current = sigma_ys_new
       Sinv_current = Sinv_new
       detval_current = detval_new
       detval_S_current = detval_S_new
       n0T=n1T
    else
       !;REJECT
       if(burn_in_done) then
         reject_sigma_ys=reject_sigma_ys+ 1
         reject_sigma_ys_vector=reject_sigma_ys_vector+ 1
       endif
    endif

  end subroutine compute_sigma_ys

  subroutine compute_tau(stepsize_tau,tau_current,reject_tau,tau_hyperparam1,&
    tau_hyperparam2,tau_pdf,dim2,y_error_t,sigma_y_current,T_indices,&
    nmeasuretotal,autocorr_vec,aa,bb,numsites,nmeasuremax,&
    detval_S_current,p0_tau,detval_current,n0T,Tinv_current,&
    detval_T_current,n0,Sinv_current,deltatime,T_new,Tinv_new,dum3,&
    C,burn_in_done)

    ! input only
    logical, intent(in) :: burn_in_done
    real, intent(in) :: stepsize_tau,tau_hyperparam1,tau_hyperparam2,aa,bb,&
      detval_S_current
    integer, intent(in) :: tau_pdf,dim2,nmeasuretotal,numsites,nmeasuremax
    real, dimension(:), intent(in) :: sigma_y_current,n0
    integer, dimension(:,:), intent(in) :: T_indices
    real, dimension(:,:), intent(in) :: deltatime,Sinv_current

    ! inout
    real, intent(inout) :: tau_current,p0_tau,detval_current,n0T,&
      detval_T_current
    integer, intent(inout) :: reject_tau

    ! out only
    real, dimension(:,:), intent(inout) :: Tinv_current

    ! really locals (temporary values)
    real, dimension(:), intent(inout) :: y_error_t,autocorr_vec,C
    real, dimension(:,:), intent(inout) :: T_new, Tinv_new, dum3

    real :: dtau, tau_new, p1_tau,n1T,detval_T_new,detval_new,pT,&
      randomu
    integer :: ii

    ! Generate new value of tau
    dtau = random_normal()
    dtau = dtau*stepsize_tau
    tau_new = tau_current + dtau

    if(tau_new .le. 0) then
       ! Skip if it generates a tau less than or equal to 0 or matrix will blow-up
       if(burn_in_done) reject_tau=reject_tau+ 1
    else

       ! Compute P1 for new value of tau
       call pdf_calc(tau_new, tau_hyperparam1, tau_hyperparam2, tau_pdf, p1_tau)

       ! create new T matrix
       do ii=1,dim2
          y_error_t(T_indices(:,ii)) = sigma_y_current(ii)
       enddo

       !$OMP PARALLEL DO default(none), &
       !$OMP             private(autocorr_vec), &
       !$OMP             shared(T_new,y_error_t,nmeasuretotal,deltatime,tau_new)
       do ii=1,nmeasuretotal
          autocorr_vec = y_error_t(ii)*y_error_t*exp(-1*(deltatime(:,ii)/tau_new))
          T_new(ii,:) = autocorr_vec
          T_new(:,ii) = autocorr_vec
       enddo
       !$OMP END PARALLEL DO

       ! compute new inverse
       Tinv_new = Ainv(T_new)

       ! calculate new value of n1T = (y-Hx)^T * Rinv * (y-Hx)

       call sgemm ('N', 'N', nmeasuretotal, numsites, nmeasuretotal, aa, Tinv_new, nmeasuretotal, &
            n0, nmeasuretotal, bb, dum3, nmeasuretotal)
       call sgemm ('N', 'N', nmeasuretotal, numsites, numsites, aa, dum3, nmeasuretotal, Sinv_current, &
            numsites, bb, C, nmeasuretotal)

       n1T = dot_product(n0,C)

       ! calculate new determinant

       detval_T_new = det(T_new)*numsites
       detval_new = detval_S_current + detval_T_new

       ! compute P1/P0
       pT=alog(p1_tau/p0_tau) - detval_new + detval_current - 0.5*(n1T - n0T)

       call random_number(randomu)              ! Generates uniformly distributed random number

       if(alog(randomu) .le. pT) then
          !;ACCEPT
          p0_tau=p1_tau
          tau_current = tau_new
          Tinv_current = Tinv_new
          detval_T_current = detval_T_new
          detval_current = detval_new
          n0T=n1T
       else
          !;REJECT
          if(burn_in_done) reject_tau=reject_tau+ 1
       endif
    endif

  end subroutine compute_tau

  subroutine compute_nu_rho(stepsize_nu,stepsize_rho,nu_current,rho_current,&
    reject_nu,nu_hyperparam1,nu_hyperparam2,nu_pdf,rho_hyperparam1,rho_hyperparam2,rho_pdf,&
    arg,distance_double,sigma_ys_current,S_new,distance,numsites,Sinv_new,nmeasuretotal,&
    nmeasuremax,aa,bb,dum3,C,n0,n0T,p0_nu,p0_rho,detval_T_current,detval_current,&
    detval_S_current,Tinv_current,Sinv_current,burn_in_done)

    ! input only
    integer, intent(in) :: nu_pdf,rho_pdf,numsites,nmeasuretotal,nmeasuremax
    logical, intent(in) :: burn_in_done
    real, intent(in)    :: stepsize_nu,stepsize_rho,nu_hyperparam1,nu_hyperparam2,&
      rho_hyperparam1,rho_hyperparam2,aa,bb,detval_T_current
    real, dimension(:), intent(in)    :: sigma_ys_current,n0
    real, dimension(:,:), intent(in)  :: distance, Tinv_current
    real(kind=8), dimension(:,:), intent(in)  :: distance_double

    ! inout
    integer, intent(inout)            :: reject_nu
    real, intent(inout)               :: nu_current,rho_current,n0T,p0_nu,&
      p0_rho,detval_current

    ! out
    real, dimension(:,:), intent(inout) :: Sinv_current
    real, intent(out) :: detval_S_current

    ! really locals (temporary values)
    real, dimension(:,:), intent(inout) :: S_new, Sinv_new, dum3
    real, dimension(:), intent(inout) :: C
    real (kind = 8), dimension(:,:), intent(inout) :: arg


    integer :: ii, jj
    real ::dnu,drho,nu_new,rho_new,p1_nu,p1_rho,ga,n1T,detval_S_new,detval_new,&
      pT,randomu
    double precision :: nu_double, rho_double,alpha
    integer (kind=4) :: nb,ize,ncalc
    real (kind = 8) :: arg_dum, k_arg_dum,ga_double

    ! Generate new value of nu and rho
    dnu = random_normal()
    dnu = dnu*stepsize_nu
    nu_new = nu_current + dnu

    drho = random_normal()
    drho = drho*stepsize_rho
    rho_new = rho_current + drho

    if(nu_new .le. 0 .OR. rho_new .le. 0 ) then
       ! Skip if it generates a nu or rho less than or equal to 0 or matrix will blow-up
       if(burn_in_done) reject_nu=reject_nu+ 1
    else

       ! Compute P1 for new value of rho and nu
       call pdf_calc(nu_new, nu_hyperparam1, nu_hyperparam2, nu_pdf, p1_nu)
       call pdf_calc(rho_new, rho_hyperparam1, rho_hyperparam2, rho_pdf, p1_rho)

       ! Create new S matrix

       nu_double = nu_new
       rho_double = rho_new
       nb = int(nu_double)
       alpha = nu_double - DBLE(nb)
       nb = nb+1

       arg = sqrt(2*nu_double)*distance_double/rho_double

       ize = 2

       call gamma(nu_double, ga_double)
       ga = ga_double

       !$OMP PARALLEL DO default(none), &
       !$OMP             collapse(2), &
       !$OMP             private(arg_dum, k_arg_dum,ncalc), &
       !$OMP             shared(S_new,sigma_ys_current,alpha,nb,ize,nu_new,ga), &
       !$OMP             shared(distance,rho_new,arg,numsites)
       do ii = 1,numsites
          do jj = 1,numsites
             arg_dum = arg(ii,jj)
             if (arg_dum==0) then
                S_new(ii,jj) = 1*sigma_ys_current(ii)**2
             else
                call rkbesl(arg_dum, alpha, nb, ize, k_arg_dum, ncalc)
                k_arg_dum = k_arg_dum*exp(-1*arg(ii,jj))
                S_new(ii,jj) = sigma_ys_current(ii)*sigma_ys_current(jj)*(1/(2**(nu_new-1)*ga)* &
                     (sqrt(2*nu_new)*distance(ii,jj)/rho_new)**nu_new*k_arg_dum)
             endif
          enddo
       enddo
       !$OMP END PARALLEL DO

       ! compute new inverse

       Sinv_new = Ainv(S_new)

       ! calculate new value of n1T = (y-Hx)^T * Rinv * (y-Hx)

       call sgemm ('N', 'N', nmeasuretotal, numsites, nmeasuretotal, aa, Tinv_current, nmeasuretotal, &
            n0, nmeasuretotal, bb, dum3, nmeasuretotal)
       call sgemm ('N', 'N', nmeasuretotal, numsites, numsites, aa, dum3, nmeasuretotal, Sinv_new, &
            numsites, bb, C, nmeasuretotal)

       n1T = dot_product(n0,C)

       ! calculate new determinant

       detval_S_new = det(S_new)*nmeasuretotal
       detval_new = detval_T_current + detval_S_new

       ! compute P1/P0
       pT=alog(p1_nu*p1_rho/p0_nu/p0_rho) - detval_new + detval_current - 0.5*(n1T - n0T)

       call random_number(randomu)              ! Generates uniformly distributed random number

       if(alog(randomu) .le. pT) then
          !;ACCEPT
          p0_nu=p1_nu
          p0_rho=p1_rho
          nu_current = nu_new
          rho_current = rho_new
          Sinv_current = Sinv_new
          detval_S_current = detval_S_new
          detval_current = detval_new
          n0T=n1T
       else
          !;REJECT
          if(burn_in_done) reject_nu=reject_nu+ 1
       endif
    endif

  end subroutine compute_nu_rho

  subroutine compute_y (yi,y_current,stepsize_y,m1,z,nmeasuremax,&
                       dum2,Dinv,n1,n0,dumx,dumy,Rinv_current,n0T,m0T,m1T,reject_y,&
                       reject_y_vector,timeindex,burn_in_done)

    ! input only
    integer, intent(in) :: yi,nmeasuremax,timeindex
    logical, intent(in) :: burn_in_done
    real, intent(in)    :: stepsize_y
    real, dimension(:), intent(in)    :: z
    real, dimension(:,:), intent(in)  :: Dinv,Rinv_current
    ! inout
    integer, intent(inout)            :: reject_y,reject_y_vector
    real, intent(inout)               :: n0T,m0T,m1T
    real, dimension(:), intent(inout) :: y_current,n0

    ! really locals (temporary values)
    real, dimension(:), intent(inout) :: m1,dum2,dumx,dumy,n1

    real :: dymod, dum,dumb,dumc,dumd,n1T,pT,randomu,y_new,m1old,dum2old
    integer :: ii

    real, external :: sdot

    ! Generate new value of y
    dymod = random_normal()*stepsize_y
    y_new = y_current(yi) + dymod

    ! subsample at observation times
    ! so only update if yi exists in the index array
    if (timeindex .ne. 0) then
      ! calculate (z-y)^T*Dinv*(z-y)
      m1old=m1(timeindex)
      m1(timeindex) = m1(timeindex) + y_current(yi) - y_new
      dum2old=dum2(timeindex)
      dum2(timeindex) = m1(timeindex)*Dinv(timeindex,timeindex)*m1(timeindex)
      m1T = sum(dum2)
    endif

    ! compute (Hx-y)^T*Rinv*(Hx-y) based on current value of y
    !
    n1(yi) = n0(yi) - dymod

    call scopy (nmeasuremax, Rinv_current(1,yi), 1, dumx, 1)
    call scopy (nmeasuremax, dumx, 1, dumy, 1)
    call sscal (nmeasuremax, n0(yi), dumx, 1)
    call sscal (nmeasuremax, n1(yi), dumy, 1)

    dum  = sdot(nmeasuremax,n0,1,dumx,1)
    dumb = sdot(nmeasuremax,n1,1,dumy,1)
    dumc = sdot(nmeasuremax,n0,1,Rinv_current(1,yi),1)
    dumd = sdot(nmeasuremax,n1,1,Rinv_current(1,yi),1)
    !Further speedup with slightly reduced numerical accuracy
    !Can replace array n1 with a scalar, only using the yi'th element
    !dumc = sdot(nmeasuremax,n0,1,Rinv_current(1,yi),1)
    !dum  = n0(yi) * dumc
    !dumd = dumc + (DBLE(n1) - DBLE(n0(yi))) * DBLE(Rinv_current(yi,yi))
    !dumb = dumd * n1

    n1T = n0T - dum + dumb - n0(yi)*dumc + n1(yi)*dumd - &
         n1(yi)*n1(yi)*Rinv_current(yi,yi) + n0(yi)*n0(yi)*Rinv_current(yi,yi)

    ! Compute P1/P0
    pT= -0.5*(m1T - m0T) -0.5*(n1T - n0T)

    call random_number(randomu)              ! Generates uniformly distributed random number

    if(alog(randomu) .le. pT) then
       !;ACCEPT
       y_current(yi) = y_new
       n0(yi)=n1(yi)
       n0T=n1T
       m0T = m1T
    else
       !;REJECT
       n1(yi) = n0(yi)
       if(burn_in_done) then
          reject_y=reject_y+ 1
          reject_y_vector=reject_y_vector+ 1
       endif
       if (timeindex .ne. 0) then
          m1(timeindex) = m1old
          dum2(timeindex) = dum2old
          m1T = sum(dum2)
       endif
    endif

  end subroutine compute_y

  subroutine compute_ratios
    use MCMC_data

    accept_vector = 1. - real(reject_vector)/real(nIt)

    accept_sigma_y_vector = 1. - real(reject_sigma_y_vector)/real(nIt)

    accept_sigma_ys_vector = 1. - real(reject_sigma_ys_vector)/real(nIt)

    accept_tau = 1. - real(reject_tau)/real(nIt)

    accept_nu = 1. - real(reject_nu)/real(nIt)

    accept_y_vector = 1. - real(reject_y_vector)/real(nIt)

  end subroutine compute_ratios

End module MCMC_subroutines

